clear all;
clc;

load mnist0v1

trainsize = 100;

X = double(data);
Y = double(label(1:trainsize));
n = size(X,1);

k = 4;

[IDX,D] = knnsearch(X,X,'K',k);

A = sparse(reshape(repmat((1:n)',1,k),[],1),reshape(IDX,[],1),reshape(D,[],1));

A = A + A';
A = (A~=0);

[S, C] = graphconncomp(A,'Directed',false);

if(S > 1) 
    error('Graph not connected! Try more neighbors.');
end

B = adj2inc(A);

L = double(B)'*double(B);
DD = diag(L);
L = 2*sparse(1:n,1:n,DD) - L;

s = 10;

[V,Deig] = eigs(L,s,'SM');
V = V(:,1:s-1);

P = repmat(Y,1,s-1).*V(1:trainsize,:);

D2 = D.^2;
IDXvec = reshape(IDX(:,2:k),[],1);
Dvec = reshape(D2(:,2:k),[],1);
Ivec = reshape(repmat((1:n)',1,k-1),[],1);

Dvec = Dvec/max(Dvec);

Q = V(Ivec,:) - V(IDXvec,:);

cvx_begin sdp
cvx_solver sedumi
    variable M(s-1,s-1) symmetric
    
    minimize (-trace(M))
    subject to
        sum(Q*M.*Q,2) <= Dvec
        M > 0
cvx_end

cvx_begin sdp
cvx_solver sedumi
    variable t
    variable gam(s-1,1)
    minimize (t)
    subject to
        P*gam >= ones(trainsize,1)
        [M gam; gam' t] >= 0
cvx_end

Ypred = int8(sign(V(trainsize+1:n,:)*gam));
calacc(Ypred,label(trainsize+1:n))